g; Feynman-Kac path integral approach 

(N 

g: for the energy spectrum of many 

1—5 

boson systems 

O ■ S. Datta 

(U ■ 

S. N. Bose National Centre for Basic Sciences, India 
^ ; February 6, 2008 

B 

' Abstract 

O , We study the ground and excited states of weakly interacing Bose 

gases ( with positive and negative scattering lengths) in connection 
with Bose Einstein condensation to test the validity of using the mean 
field theory and Born approximation. They behave as new quantum 
fluids ( a gas in the weak limit and a liquid in the dense limit and we 
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I study their many body physics in the dilute limit within the realistic 

■ potential model ( Morse type) by Feynman-Kac path integral tech- 

' nique . Within numerical limitations, this method is exact in priciple 
and turns out to be a better alternative to GP as all the ground and 

^ ' excited states properties can be calculated in a much simpler way. 
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1 Introduction 



With the experimental reahzation of Bose Einstein Condensation in alkah 
gases[l], the study of many boson systems has become an area of active re- 
search interst. Previous numerical procedures are based on mean field theory 
hke Gross Pitaevski [2] etc. They seem to work well for ground state prop- 
erties but turn out to be approximate as they fail to include correlations in 
the many body theory. Investigations of effects beyond mean field theory is 
an important task and makes the many boson systems interesting even from 
the many body perspective [3]. Moreover earlier calculations with S function 
potential do not solve the many body problem exacxtly but only within a 
perturbation theory as in a system of Bose gas with 6 function potential, 
particles do not collide [4]. As a result, speculations ( particularly excitation 
frequencies etc ) based on these methods differ drastically for different exper- 
imental modes. So an alternative to GP was necessary which can describe 
the effect of interaction in a more reliable way and predict the excitation 
frequencies and other properties more accurately. Eventhough Monte Carlo 
techniques are slow, computationally expensive and faces sign problem for 
fermionic systems, these are the only numerical techniques available for these 
kind of many body systems, which are exact and can include corrlations in a 
reliable way. We also test the validity of Born approximation at low energy 
nd temperature. 

We propose to apply a quantum Monte Carlo technique based on 
Feynman-Kac formalism [5-7] of Quantum Mechanics which forms the basis 
for a simple and accurate way to calculate ground and excited state prop- 
erties. To be precise, we apply Generalized Feynman Kac ( GFK) method 
[8-10] for attaractive and repulsive potentials for many boson system of al- 
kali gases ( Rb and Li respectively ).In our model the atoms in the dilute gas 
interact through Morse potential instead of conventional pseudopotential. 
Since at low temperature the de Broglie wavelength of the atoms become 
appreciable, we do full quantum treatment. To connect GFK to other many 
body techniques our numerical procedure ( GFK ) has a straightforward im- 
plementation to Schroedinger's wave mechanics. GFK is essentially a path 
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integral technique with trial functions for which operations of the group of 
the wave function keep points in the chosen nodal region, provide an up- 
per bound for the lowest state energy of that symmetry. The nodal region 
with the lowest energy serves as a least upper bound. If the nodal region 
has exact nodal structures of the true wave function the random walk is 
exact in the limit scale, time for walk, and number of walks get arbitrarily 
large. Rh^^ and Li' both having odd number of electrons and odd num- 
ber of nuclei obey Bose statistics. So we do not need to worry about sign 
problem for these calculations. Our method gives more accurate values for 
the ground and excited state energies as to calculate energy we approximate 
an exact solution ( i.e., the GFK representation of it ) to the Schroedinger's 
equation, whereas most of the other numerical procedures approximate a 
solution to an approximate Schroedinger equation. Since we work with the 
bosonic system our results are exact within the numerical limitations. There 
are other Monte Carlo calculations which are worth mentioning in this con- 
text are ground state calculation by Difiision Monte Carlo method [lljand 
path integral Monte Carlo ( with temperature dependence )[12]. This pa- 
per is organized as follows. In Sec 2, we present the general formalism of 
Feynman-Kac and Generalized Feynman-Kac method and then describe it 
for the ground and excited state of trapped Bosons. In Sec 3, we describe 
the numeraical procedure. In Sec 4 we present results for positive ( it!6^^ ) 
and negative { Li^ ) scattering lengths. Finally in Sec 5 we summerize our 
achivements. 



2.1.1 Feynman-Kac Path integretion 

For the Hamiltonian H — —A/2 -|- V{x) consider the initial value problem 
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Path integral Theory at T=0 



.du 





(1) 
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with X & R'^ and u{0,x) = 1. The solution of the above equation can be 
written in Feynman-Kac representation as 

u(t, x) = E^exp{- f V{Xis))ds} (2) 
Jo 

where X(t) is a Brownian motion trajectory and E is the average value of 
the exponential term with respect to these trajectories. The lowest energy 
eigenvalue for a given symmetry can be obtained from the large deviation 
priniciple of Donsker and Varadhan [13], 

A = - lim -lnE^exp{- f ViXis))ds] (3) 

t-»oo t Jo 

Generahzations of the class of potential functions for which Eqs. 2 and 3 
are valid are given by Simon[14] and include most physically interesting po- 
tentials, positive or negative, including, in particular, potentials with 1/x 
singularities. It also means that the functions determined by Eq(2) will be 
the one with lowest energy of all possible functions independent of symme- 
try. Restrictions on allowed Brownian motions must be imposed to get a 
solution of the desired symmetry if it is not the lowest energy solution for a 
given Hamiltonian. Although other interpretations are interesting, the sim- 
plest is that the Brownian motion distribution is just a useful mathematical 
construction which allows one to extract a physically relevant quantity, the 
ground and excited state energy of a quantum mechanical system. In nu- 
merical implementation of Eq(3) the 3N dimensional Brownian motion is 
replaced by 3N independent, properly scaled one dimensional random walks 
as follows. For a given time t and integers n and 1 define [6] the vector in 
R^^ 

W{1) = W{t,n,l) = {wi\t,n,l),W2\t,n,l),Ws\t,n,l).... (4) 

Wi'^{t, n, l)w2^ {t^ n, l)ws^{t, n, I) 

where 

wi{t,n,l) = i:^ (5) 
k=i 

with ^^'(0, n,l) — for i = 1, 2, ....N;j — 1, 2, 3 and I — 1, 2, , , , , , nt. Here 
e is chosen independently and randomly with probability P for all i,j,k such 
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that -P(ejfe — 1)—P{e]ji. — —1)—^ It is known by an invariance principle[15] 
that for every u and W(l) defined in Eq(4) 

-| nt 

il^^(-Em(0))<^ (6) 

^ 1=1 

t 

^P{j V{X{s))ds < V 



Consequently for large n, 

t 

P[exp(- j V{X{s))ds) < v] (7) 



1 nt 

^ P[exp( ^(^(0)) < ^] 

^ 1=1 

By generating iV^ep independent replications Zi,Z2,....Zn^^^ of 

1 nt 

z^^eM-i--j:nwm (s) 

1=1 

and using the law of large numbers, {Zi + Z2 + ...Z]si^^^)/Nrep — Z{t) is an 
approximation to Eq(2) 

A — logZ{t) (9) 

t 

Here W^™(/),m = 1,2, N^^p denotes the m*^ realization of W(l) out of iV^ep 
independently run simulations. In the limit of large t and A^^ep this ap- 
proximation approaches an equality, and forms the basis of a computational 
scheme for the lowest energy of a many particle system with a prescribed 
symmetry. 

In dimensions higher than 2, the trajectory X(t) escapes to infinity 
with probability 1. As a result, the important regions of the potential are 
sampled less and less frequently and the above equation converges slowly. 
Now to speed up the convergence we use Generalized Feynman-Kac (GFK) 
method. 

2.1.2 Generalized Feynman Kac path integretion 

To formulate the (GFK) method, we first rewrite the Hamiltonian as H — 
Ho + Vp, where Hq = -A/2 + Xt + AipT/^ipT and Vp = 1/- (At + A^t/2^/^t). 
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Here tpT is a twice differentiable nonnegative reference function and Hi/jt — 
XxipT- The expression for the energy can now be written as 

A = At - hm -lnE^exp{- t VJY{s))ds} (10) 

t^oo t Jo 

where Y(t) is the diffusion process which solves the stochastic differential 
equation 

The presence of both drift and diffusion terms in this expression enables the 

trajectory Y(t) to be highly localized. As a result, the important regions of 
the potential are frequently sampled and Eq (3) converges rapidly. 



2.2 Schroedinger Formalism for Rb and Li condensate 

Case 1: A system of few Rb atoms interacting through a model potential in 
an isotropic trap. 



1 ^ 

[-A/2 + V^nt + 2 E[^i' + y^'' + = (12) 

i=l 

The atoms here interact through a model potential Dsech'^{r/ro) 
Case 2 : Systems of Rb and Li atoms interacting through Morse potential 
trapped inside an anisotropic trap. We choose to work in the cylindrical 
coordinates as the original experiment had an axial symmetry, cylindrical 
coordinates are the natural choices for this problem. We consider a cloud of 
N atoms interacting through repulsive potential placed in a three dimensional 
anharmonic oscillator potential. At low energy the motion of condensate can 
be represented as 

1 ^ 

[-A/2 + V,^t + t; 121^' + y^' + A^.Xf) = Et^(f) (13) 
^ i=i 

where | Y^^i[xi'^ + yi'^ + Xzi'^]tlj{f) is the anisotropic potential with anisotropy 
factor A = ^. Now 

Vint = VMorse = E ^(^i.) = E D[e-^^'-''^^ (e-«(-'-o) - 2)] (14) 
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The above Hamiltonian is not separable in spherical polar coordinates be- 
cause of the anisotropy. In cylindrical coordinates the noninteracting part 
behaves as a system of noninteracting harmonic oscillators and can be writ- 
tem as follows : 

\- — — ( 

+i(/ + AV)]^(p,^) 
= E^l^{p,z) (15) 

The energy 'E' of the above equation can be calculated exactly which is 

En.n.m = (2/1^ + \m\ + 1) + {u, + l/2)A (16) 

Since we are dealing with many Boson systems having very small 
number of particles, ( In JILA experiment the number of particles is of the 
order of 2000 and does not correspond to Thomas- Fermi limit ). So it is quite 
legitimate to use Gaussian trial functions for modes which are not coupled. 
In our guided random walk we use the solution of Schoroedinger equation for 
harmonic oscillator in d-dimension as the trial function as follows [16]: 

^n,„.„^(f) ^ exp^ Hr^Sz) X e'^^^ p^e'^' L^}'^\p^) (17) 
2.3 Fundamentals of BEC 

Even though the phase of Rb vapors at T=0 is certainly solid, Bose con- 
densates are preferred in the gasous form over the liquids and solids because 
at those higher densities interactions are complicated and hard to deal with 
on an elementary level. They are kept metastable by maintaining a very 
low density. For alkali metals, r], the ratio of zero point energy and molec- 
ular binding energy lies between 10~^ and 10~^. According to the theory of 
corresponding states[17] since for the T=0 state of alkah metals, 77 exceeds 
a critical value 0.46, the molecular binding energy dominates over the zero 
point motion and they condense to solid phase. But again the life time of a 
gas is limited by three body recombination rate which is proportional to the 
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sqaure of the atomic density. It gets suppressed at low density. Magnetically 
trapped alkali vapors can be metastable depending on their densities and life- 
times. So keeping the density low only two body collisions are allowed as a 
result of which dilute gas approximation [11] still holds for condensates which 
tantamounts to saying na^ << 1 (a is the scattering length of s wave). Now 
defining n = N/V = as a mean distance between the atoms ( definition 
valid for any temperature ), the dilute gas condition reads as a << VaV and 
zero point energy dominates (dilute limit). In the dense limit, for a ~ Vav 
on the other hand the interatomic potential dominates [11] .The gas phase is 
accomplished by reducing the material density through evaporative cooling. 

3 Numerical procedure 
3.1 Dilute limit 

In the dilute limit and at very low energy only binary collisions are possible 
and no three body recombination is allowed. In such two body scattering at 
low energy first order Born approximation is applicable and the interaction 
strength 'D' can be related to the single tunable parameter of this problem, 
the s-wave scattering length 'a' through the relation given below. This single 
parameter can specify the interaction completely without the detail of the 
potential in the case of pseudopotentials. We study a system of Rb particles 
interacting through following two repulsive potentials. 

A. Model potential V{r) = Dsech^{r/ro) 

B. Morse potential for dimer of rubidium 



with more realistic feature of having repulsive core at = than other 
model potentials. In our case, the interaction strength depends on two more 
additional parameters, Tq and a. 



^V{rij) = ^D[e-"(''-'-o)(e-°(^-''°) - 2)] 



(18) 




(19) 
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Case I. For potential 'A', the above yields D — ^^"^^3 The Hamiltonian for 
a system of Rb gas in an isotropic trap intertacting throuh a potential V{r) 
can be written as 

N 



/2 
i 

i=l i,j 



1=1 

= E^{j^) (20) 

The above Hamiltonian can be rescaled by substituting f — sr and E — EqU 
and the interaction strengh in the dimensionless form can be written as 

12a o 
71 = — = 3.05 X 10-=^ 21 

with To = 1.2. 

Case II. Again for potential 'B' the Born approximation yields 

D = ; (22 

^gQro(-garo _ 16) ^ ' 

The Hamiltonian for Rb gas with an asymmetric trapping potential and 
Morse type mutual interaction can be written as 

N 

i=l i,j 

N N N 

E + E y'^ + E 

i=l i=l i=l 

= EiP{f^) (23) 

The above Hamiltonian can be rescaled by substituting f* = sr and E — EqU 
as 

^ 2ms2^ ^ ^^ms3e-o(e-o_i6)L^ 
ms"^ ^ ^ ^ 

E + E y^' + E 

^ i=\ i=l i=l 

= EoU^{r) (24) 
, jv 3 

I2 g - ^e°-.(e°°. - 16) Eie-'-'-'^He-"'"'-' - 2)] 
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2 2 4 iV 2 2 

m UJx S ^ , 2 2 '^z 2 



= -Eo^^P{r) (25) 

Let niW^ = 1 =^ s2 ^ ^ is the natural unit of length. Let = 1 =^ 
[/ — = ^a;^; is the natural unit of energy. Then the standard form of the 
equation becomes 

[-T^i" - T4 y[e-^(rij-ro)f^^-a{n,-ro) _ 3)] 

^ j=l '^x ^x 

= -£;oV'(f) (26) 
With ujx — ujy = the above eqn becomes, 

^2^^^ "^56-0(6-0 - 16) ^ 

1=1 V / l<] 

1 ^ 

i=l 

= -Eoi^i^ (27) 



1 ^ 

hY.^i^ - ^^[e-(---'o)(e-"('-^^-^o) - 2)] 

^ i=l i<j 
1 ^ 

^ i=i 

= -£^oV'(0 (28) 

Now for a — .35 and Tq = 11.65 (both in oscillator units)[4,18], the Morse 
potetial become almost noninteracting and the results become substantially 
lower than GP. We find a better agreement with GP by choosing a — .29 
and To = 9.67. We have checked that for these choice of parameters, Morse 
solution is extremely good, a = 52 x 10~^°, s = .12 x 10^^, the interaction 
strength 7 is given by 

72 = 4 = 4.9 X 10-^ (29) 
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For mean field calculation the value of interaction strength was taken to 
be 4.33 X 10~^. For this problem we are interested in the limit 7 << 1. 
The case 7 >> 1 is usually known as the Thomas Fermi limit. For 7 = 
4.9 X 10~^, the eigenvalue equation reduces to a minimally perturbed system 
of d dimensional anisotropic oscillator where d — 3N and N is the number 
of particles. The whole concept of bound states of Morse dimers is very 
outside the range of this limit, so the nonexistence of two-body bound states 
is ensured by choosing the above parameters. 

Case III. Using the same Morse potential with a — 0.4 and Tq = 
6.8, one can generate attractive interaction for as well. In this case the 
interaction strength turns out to be 

73 = 4 — = -9.345 X 10-6 (30) 

Even though 7 << 1, we solve the eigenvalue cqn nonpcrturbatively with 
Generalized Feynman-Kac procedure. Energies at zero temperature are ob- 
tained for ground and excited states by solving Eq.(l) and using Eq.(3). 
Since original Feynman-Kac method [5-6] is computationally inefficient we 
incorporate importance sampling in our random walk and use trial function 
of the form given in Eq.(17) 
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4 Results 



Case 1. A system of Rb atoms in an isotropic trapping and interactinng 
through a model potential 
Dsech'^{r /tq). 

Table 1: In column 2, ground state energies for 3-20 particles in an isotropic 
trap with potential 'a' in units of 1%uj±_ for a — .00305a/io- Columns 3-5 
represent the same energies by Blume and Green(BG)[19], GP and modified 
GP respectively 



N EqfK EbG Eqp Eop^mod 

3 4.50925(3) 4.51037(2) 4.51032 4.51032 

5 7.53086(5) 7.53439(6) 7.53432 7.53434 

10 15.13867(3) 15.1539(2) 15.1534 15.1535 

20 30.58460(4) 30.639(1) 30.638 30.639 



4.1 Positive scattering length: Rb 

As the potential does not sustain any many body bound state ( which is 
ensured by suitably choosing the value of the parameters a and ro ) and the 
scattering length is positive the system behaves as a gas or as a metastable 
state which can be long-lived at very low densities[3]. In the table below, we 
explicitly show the expectation values of trap potential, interatomic potential 
and kinetic energy as three components of total energy for different number 
of particles and it is observed that virial theorem is satisfied in each case[20]. 

"^Ekin — 2Eho + ^Epot — (31) 

Prom Pig 1, 2 and 3, we see that energy /particle rises with increase in number 
of atoms in the trap for different symmetry states. Prom Pig 4,5 and 6, 
we see that in the case of Rb, increase in number of particles in the trap 
lowers the central density. Prom the ground state data,the aspect ratio of 
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velocity distribution i.e., y^^f^ is calculated and its values are given in 
the last column of Table 2. The value of aspect ratio in the noninter acting 
limit corresponds to -\/A and with increase in number of particles this value 
increases and in the Thomas Fermi limit it should correspond to A. Prom 

our data, we verify this trend of aspect ratio( Fig 7 ). 

Table 2: Results for ground state of Rb with A = ^/S Chemical potential 
and energy are in units of huj± and length is in units of a±. Numbers in the 
brackets correspond to the reference[20] 



N /i E/N E/N,^^ E/N,^ E/N^^^ ^ < > / < > 



1 


2.414213 


2.414213 


1.207409 


1.206803 


0.0 


1.679563 














(1.68333) 


10 


2.448952 


2.431595(5) 


1.202488 


1.211725 


.017369 


1.684900 


40 


2.564350 


2.489287(1) 


1.196455 


1.217758 


0.075068 


1.68732 


70 


2.678893 


2.546602(6) 


1.188591 


1.225621 


0.132339 


1.688079 


100 


2.792482 


2.603549(3) 


1.180656 


1.233556 


0.189134 


1.688960 




(2.88) 


(2.66) 








(1.79545) 


200 


3.149535 


2.79075(7) 


1.061734 


1.35247 


0.367660 


1.690056 




(3.21) 


(2.86) 








(1.88888) 
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Figure 1: A plot for the Condensate Energy/ Particle versus Number of atoms 
in trap for 200 particles for the ground state : Uz — rip — m — 0; this work 



Table 3: Results for excited states of Rb gas with A = VS 



N 


/ nz=np=0,m=l 


^ 1 nz=np=0,m=2 


^ / nz=np=0,m=3 


10 


3.428526(1) 


4.424937(2) 


5.422566 (1) 


40 


3.474993(3) 


4.462418(3) 


5.445864(4) 


70 


3.520849(5) 


4.500300 (4) 


5.472183(6) 


100 


3.566892(4) 


4.537793 (6) 


5.494076(7) 


200 


3.718571(8) 


4.662026(7) 


5.590799(8) 
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50 100 150 200 

Number of Particles [ N ] 

Figure 2: A plot for the Condensate Energy/Particle versus Number of atoms 
in trap for 200 particles for the 1st excited state : Uz — rip — 0,m — l;this 
work 




4.4 I , , \ , , \ , , \ 

50 100 150 

Number of Particles [ N ] 



Figure 3: A plot for the Condensate Energy/ Particle versus Number of atoms 
in trap for 200 particles for the 2nd excited state: Uz — Up = 0,m = 2; this 
work 
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Figure 4: Ground state wave function of Rb along the x axis for N = 10 [this 
work] 




Figure 5: Ground state wave function of Rb along the x axis for N = 40 [this 
work] 
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Figure 6: Ground state wave function of Rb along the x axis for N = 100 
[this work] 



03 
O 



50 



100 

N 



150 



200 



Figure 7: Aspect ratio in Rb^7 as a function of N. The horizontal line corre- 
sponds to a number close to \/X [this work] 
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4.2 Negative scattering length: Li 

In this section we report the calculations for Li. Since Li has negative scat- 
tering length, in Li gas the atoms interact with each other with attractive 
force and when this force becomes larger than molar binding energies the 
gas collapses. But it is observed that if the number of particles in the gas 
is not too high, then zero point energy still can exceed this attaractive force 
and still there can be a metastable state. Since we are dealing with only 
200 atoms this does not pose any threat to our work. As Morse potential 
can be used both for positive and negative scattering lengths [21], we a use 
Morse potential for Li gas too and choose a = 0.4 and ro = 6.8 to produce 
an attractive type of interaction between Li atoms. In this case anisotropy 
factor is A = 0.72 and scattering length is — 27ao- From Fig 8, 9 and 10, it 
is evident that increase in number of particles sharply increases the central 
density with contrast to Rb gas with positive scattering length where it de- 
creases with increase in number of particles. In Fig 11, we plot aspect ratio 
of Li gas as a function of number of particles in the gas. 
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Table 4: Results for ground state of Li, A = -s/S Chemical potential and 
energy are in units of hco± and length is in units of a±. Numbers in the 

brackets correspond to the reference [22] 



N 


A* 


E/N 


' Kin 


1 no 


' poz 


J< x'^ > / < z'^ > 
V ' 


1 


1.360000 


1.360000 


0.680273 


0.679726 


0.00000 


0.847555 






1.360000 








(0.84873) 


10 


1.354930 


1.357467(1) 


0.680486 


0.6795132 


-0.002534 


0.848088 


40 


1.338018 


1.349014(3) 


0.681689 


0.678310 


-0.010990 


0.848529 


70 


1.32108 


1.340545(5) 


0.682509 


0.677490 


-0.019458 


0.848998 


100 


1.304113 


1.332061(6) 


0.683446 


0.676553 


-0.027943 


0.849123 




(1.327) 


(1.344) 








(0.850728) 


200 


1.247306 


1.303674(7) 


0.687241 


0.672758 


-0.056346 


0.849731 




(1.291) 


(1.326) 








(0.854858) 



^ 1.2245 




Figure 8: Ground state wave function of Li along the x axis for N = 10 [this 

work] 
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Figure 9: Ground state wave function of Li along the x axis for N = 40 [this 

work] 



^ 5.5e+09 - 



0.44 
X 



Figure 10: Ground state wave function of Li along the x axis for N = 100 
[this work] 
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1.8475 ' ' ' ' ' ' ' ' ' ' ' ' ' 

50 100 150 200 

Number of atoms [ N ] 

Figure 11: Aspect ratio in If as a function of N. 



5 Conclusions: 

Numerical work with bare Feynman-Kac procedure employing modern com- 
puters was reported [6] for the first time for few electron systems after forty 
years of original work [5] and seemed to be real useful for calculating atomic 
ground and excited states [7]. Tremendous success [9] in atomic physics mo- 
tivated us to apply it to condensed matter physics. We have performed our 
calculations only with 200 particles. Nonetheless we have been able to demon- 
strate some of the 'Holy Grails' of Bose Condensation, namely we have been 
able to calculate energy, chemical potential, aspect ratio etc which reason- 
ably agree with existing results achieved through other methods. The reason 
for this is that even with small number of particles we have sovled many 
boson systems nonpcrturbatively and fully quantum mechanically. Within 
the numerical limitations our many body technique is exact as solve the full 
Hamiltonian with reahstic potentals. We have calculated spectrum of Rb 
and Li gas by considering realistic potentials like Morse potential, instead of 
conventional pseudopotentials for the first time, thus allowing one to do the 
calculation exactly as opposed to the case of S function potential where it is 
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calculated perturbatively. 

We have been successful in achieving a lower value for Rb in both 
isotropic and anisotropic cases and Li ground state than that obtained by 
Gross-Pitaevski technique [2] . As a matter of fact in all three cases our re- 
sults stay lower than those coming from GP. For the isotropic case ( Case 
1 with a system of 3-20 Rb atoms ) this does not contradict results of Lieb 
and Yngvason (LY) [23] which states that energies obtained from GP equa- 
tion provide an lower bound to the actual ground state of a Bose gas with 
non-negative, finite range, spherical two-body potential. This is due to the 
fact that mean field theory inherent in GP equation does not apply for sys- 
tem with small number of particles [11]. In the case of anisotropic case, LY 
theorem does not apply as theinteracting two body potential is non spherical 
in nature. Only general principle that applies here is the energies coming 
from the GP technique used in ref [1] are variational estimates (results in this 
case are obtained by minimizing the energy functional corresponding to the 
Hamiltonian with mean field interaction instead of having some correlations 
explicitly.) and provide an upper bound for the diffusion Monte Carlo or 
similar exact calculations. This is again in conformity with Dyson's upper 
bound [24] for low density Bose gas. We also found that Born approximation 
is valid at low energy and low temperature. Using this approximation we 
could produce results in all three cases which favorably compare to those 
given in the litearture. However it is crucial to choose the paramters used in 
the potential correctly as these parameters control the strength of the inter- 
action directly. There are no hard and fast rules for how weak the intercation 
should be as it was pointed out in the noble lecture[25]. " It is far more di- 
lute and weakly interacting than liquid helium super fluids, for example, but 
far more strongly intercating than the nonintercating light in a laser beam" . 

We have found an alternative to Gross-Pitaevski technique and 
other mean field calculations. The method is extremely easy to implement 
compared to mean field tecniques and our fortran code at this point consists 
of about 270 lines. The simplicity in our many body theory is appealing as 
the mere ability to add, subtract and toss a coin enables one to solve the many 
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body problem. We employ an algorithm which is essentially parallel in nature 
so that eventually we can parallelize our code and calculate thermodynamic 
properties of bigger systems ( of the order of 2000 atoms ) taking advantage 
of new computer architechtures. This work is in progress. We are continuing 
on this problem and hope that this technique will inspire others to do similar 
calculations. 
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